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SOME EXPERIENCES WITH THE 
VISCOUS-INVISCID INTERACTION APPROACH 


W. R. Van Dalsem, J. L. Steger 
NASA Ames Research Center, Moffett Field, California 

and 

K. V. Rao 

Iowa State University, Ames, Iowa 

Methods for simulating compressible viscous flow using the viscous-inviscid interac- 
tion approach axe described. The formulations presented range from the more familiar 
full-potential/boundary-layer interaction schemes to a method for coupling Euler/Navier- 
Stokes and boundary-layer algorithms. An effort is made to describe the advantages and 
disadvantages of each formulation. Sample results tire presented which illustrate the ap- 
plicability of the methods. 


NOMENCLATURE 
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Cf 

Cl 

C, 

c 

Cp > 

E 

H 

h 

M 

Pr 

P 

R 

Re 

s 

T 

t 

U,V,W 

UyVyW 

*,y,z 


vector potential 
skin friction coefficient 
lift coefficient 
pressure coefficient 
chord 

specific heats at constant pressure and volume, respectively 

shift operator (e.g., E^Uj = Uj± i) 

total enthalpy or shape factor 

enthalpy or time step 

Mach number 

Prandtl number 

pressure 

perfect gas constant 
Reynolds number 
entropy 
temperature 
time 

contravariant velocity components 
velocity components 
coordinate axes 


0 (P/PUD 

7 ratio of specific heats 

A forward difference operator (e.g., A^ = (E^ 1 — 1 )/(A£)) 

V gradient or backward difference (e.g., V( = (1 — E^' 1 )/(A^)) operator 

8 central difference operator (e.g., 8{ = (E^ 1 — E^" 1 )/( 2 A^)) 
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6 

6 
l >* 

0 

6 

K 

P 

t,V,( 

P 

<t> 

u> 

<jj 


upwind difference operator, see §B.2 of Part I 

midpoint central difference operator (e.g., l )/(A«) 

displacement thickness 

flow deflection angle 

momentum thickness 

total conductivity 

total viscosity 

transformed coordinates 

density 

scalar potential function 
vector potential functions 
vorticity vector 

overrelaxation parameter or vorticity magnitude 


Superscripts /Subscripts 


e 

edge of boundary layer 

i 

inviscid value 


£, 77 , and £ indices 

71 

iteration number 

r 

reference value 

t 

derivative with respect to time 

V 

viscous value 

w 

wall value 

0 

stagnation value 

oo 

free stream value 


INTRODUCTION 

The simulation of three-dimensional flows by solving finite-difference approximations 
to the Navier-Stokes equations is becoming more routine. However, because of large com- 
puter memory and time requirements, it is often impractical to compute these flows using 
the fine grids necessary to resolve all of the high gradients and accurately predict some 
global parameters, such as the drag. Improving accuracy by grid refinement can be costly 
because more iterations are required to converge to a solution and each iteration costs 
more. Increasing the resolution by an order of magnitude may, for example, increase the 
required computer time by two orders of magnitude. 

In many applications the Navier-Stokes equations are solved throughout the entire flow 
field, although the majority of the flow can be accurately described by simpler equations 
which can be solved using less computer resources. An alternative is to take advantage of 
the physics of the specific zones of the flow. For example, the potential flow far from the 
body can be resolved by solving the scalar potential equation, shear layers can be simulated 
by solving the boundary-layer equations, while the Euler or Navier-Stokes equations can 
be used in any remaining regions. The savings will be significant if the optimal algorithm 
is used to accurately resolve each flow region. However, a potential disadvantage of such a 
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scheme is the increased complexity (and possible numerical instability) of treating multiple 
solution algorithms, solution domains, and domain interfaces. This difficulty is usually 
accentuated when the generality of the Navier-Stokes equations must be retained. 

This paper describes some of our experiences in applying viscous-inviscid interaction 
methods as well as efforts towards generalizing and simplifying these procedures. Because 
efficient boundary-layer algorithms are the “backbone” of all viscous-inviscid interaction ef- 
forts, two- and three-dimensional finite-difference boundary-layer algorithms are described 
in Part I of this paper. Some of the viscous-inviscid interaction schemes which we have 
explored are described in Part II. 

PART I: BOUNDARY-LAYER FORMULATIONS 

Many important flows include mildly separated boundary layers which still satisfy the 
thin shear-layer assumption used in deriving the boundary-layer equations . 1 Therefore, it 
was natural that early investigators tried to use the boundary-layer equations to study 
separation. However, it was soon realized 3 that any attempt to march the boundary- 
layer equations through a separation point met with failure. This motivated Goldstein 3 to 
study the behavior of the boundary-layer equations at separation. Goldstein showed that 
near separation the boundary-layer equations may generate a singular solution unless the 
prescribed pressure distribution satisfies very restrictive relations. Hence, researchers spec- 
ulated that the pressure distribution should not be merely specified to the boundary-layer 
equations, but should result from a strong interaction between the viscous and inviscid flow 
fields . 4 This concept resulted in the successful application of the boundary-layer equations 
to separated flow. 

Two approaches arose for marching the boundary-layer equations through a sepa- 
ration point. The first to appear was based on the concept of solving the inviscid and 
boundary-layer equations simultaneously. This approach was first applied to the inter- 
action between a supersonic inviscid flow and a separating boundary layer . 5-9 The con- 
cept has since been extended to the computation of the transonic separated flow about 
airfoils , 10 and the computation of steady and unsteady diffuser flows , 11-13 to cite two ex- 
amples. A variation of this concept known as the quasi-simultaneous approach has also 
been used. In this approach the full inviscid equations are solved, and the viscous flow is 
then computed simultaneously with a simplified inviscid flow model. This cycle continues 
until convergence . 13 However, it has proven challenging to extend either the simultaneous 
or quasi-simultaneous viscous-inviscid interaction concepts to allow the use of the more 
complex inviscid and viscous flow models in external flow applications. 

The second method of applying the boundary-layer equations to separated flow arose 
from the observation that if the pressure was not specified to the boundary-layer equations, 
nonsingular separation point solutions would be obtained. Of course, some forcing function 
other than the pressure distribution would have to be specified. Catherall and Mangier 14 
successfully marched the boundary-layer equations through a separation point by specify- 
ing a displacement thickness rather than a pressure distribution. Since the displacement 
thickness is typically one of the results of a boundary-layer computation, Catherall and 
Mangler’s algorithm was referred to as an inverse method, in contrast to the classical 
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direct pressure specified algorithms. Both integral 15 16 and finite-difference 17-21 inverse 
boundary-layer algorithms have since been developed. 

Solution algorithms for the boundary-layer equations are perhaps equally divided be- 
tween those which use integral methods and those using either finite-difference or finite- 
element methods. Integral boundary-layer methods are quite fast and with use of good 
correlation functions, they can give good accuracy. However, they require modeling of 
the flow physics in addition to the turbulence model. As a result it is difficult to de- 
cide whether errors axe due to turbulence modeling or to the integral method modeling. 
The integral methods are also user intensive since they require extensive correlation with 
experimental data. Finally, accurate integral methods have not appeared for separated 
three-dimensional flow. For these reasons we have preferred finite-difference boundary- 
layer methods. (It has also been our experience that once the boundary-layer and inviscid 
flow equations axe coupled, the computer time savings accrued by using an integral over a 
finite-difference boundary-layer scheme is minimal.) 

A. Steady Two-Dimensional Formulation 

1 . Equation Set 

The nondimensionalized, compressible boundary-layer equations for the steady, two- 
dimensional flow of a perfect gas are 
x-momentum equation 

p(UU X + VUy) = -(3p x + ( flUy) v (id) 

energy equation 

pc p (uT x + vT v ) = fiup x + (/cT y ) y + p(u v ) 2 (lb) 

perfect gas equation 

P = P T ( lc ) 

continuity equation 

+ (pv)v = 0 ( ld ) 

where x,y are along and normal to the body or wake centerline, respectively, and u and v 
are velocity components in the corresponding directions. Viscosity, pressure, temperature, 
and density are nondimensionalized by their free stream values. The u and x variables axe 
nondimensionalized by the free stream velocity and a characteristic length, respectively. 
The v and y variables are nondimensionalized by the same quantities divided by y/Re^. 

To solve the boundary-layer equations using finite-difference approximations it is nec- 
essary to construct a grid. In the interest of accuracy and computational efficiency it is 
desirable to use a nonuniform grid. However, if Eqs. (la-d) are differenced on a nonuniform 
grid, complicated variable-spacing finite-difference operators must be used. This problem 
can be avoided if a general x,y to £(x), 77 (x,y) coordinate transformation is applied to the 
boundary-layer equations. A nonuniform x,y grid may then be used while the equations 
can be solved on a uniform £,77 grid. Thus, the physical domain x,y grid is adapted to 
resolve the flow field, while the computational domain £,77 grid is chosen so that simple 
equal-spaced finite-difference operators may be used (Fig. 1). Although the equations are 
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slightly more complicated after the coordinate transformation, overall the finite-difference 
solution of the equations is much simpler. 

For computational convenience, Eqs. (la-d) are transformed from the physical x,y 
domain to a uniform £(z),t/(x,t/) computational domain 
x-momentum equation 

/>[«(«*& + UvVz) + vu n rty] = -PptZ* + (2 a) 

energy equation 

pCp[u(T(£ x + T v r] z ) + vT^rjy] — + ( K Ti)T}y)r,T]y + fx(u n rj y ) 2 (2b) 

perfect gas equation 

P — pT (2c) 

continuity equation 

(/>«)(£* + (pu)vVz + (pv) v Vv = 0 (2d) 




COMPUTATIONAL 

DOMAIN 
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Fig. 1 Two-dimensional boundary-layer physical to computational domain co- 
ordinate transformation. 


2. Solution Algorithm 

Because the steady boundary-layer equations are parabolic in space, it is convenient 
and economical to solve the equations by marching along the body. In two-dimensional 


5 



flows, the marching direction is clear and space marching schemes are quite successful. 
However, in three dimensions the marching direction is not clear and there are difficulties 
in applying marching schemes to complex three-dimensional flows. (We will return to this 
point in the next section.) 


The relevant length scales in a boundary layer (e.g., the boundary-layer thickness and 
a characteristic length in the flow direction) differ in order by a power of the Reynolds 
number. For high Reynolds numbers, this implies a stiff numerical problem and implicit 
schemes are generally required. One of the more popular of these algorithms is the “box 
scheme” originally developed by Keller and Cebeci, 22 which Cebeci has shown is appli- 
cable to a wide range of thin shear-layer problems. 23 However, this and other implicit 
algorithms which solve the boundary-layer equations simultaneously require the solution 
of block tri diagonal matrices. Since the boundary-layer equations are weakly coupled, a 
numerical algorithm which solves the boundary-layer equations in a sequential manner at 
each stream wise station may be more appropriate and efficient. 


The boundary-layer equations can be solved using an implicit predictor-corrector 
algorithm. 24-25 Streamwise marching begins by predicting estimates of u, T, p, v, p, and 
k . This predictor step uses first-order-accurate £ operators in the x-momentum and energy 
equations, but for one step produces second-order-accurate values. Because the predictor 
step only needs to be first-order accurate, any nonlinear coefficients can be lagged in 
£. A second-order-accurate corrector step is then used to calculate improved values of 
u, T, p, u, p, and k . During the corrector step, the nonlinear coefficients are evaluated 
using the most recently computed flow variables. At each step the equations are solved 
sequentially. First u is updated from the x-momentum equation. With u available, T is 
updated from the energy equation, p is then obtained from the equation of state, and v 
can then be updated from the continuity equation. Overall, with this predictor-corrector 
scheme, second-order-accurate solutions are obtained at the cost of two scalar bidiagonal 
and four scalar tridiagonal matrix inversions per streamwise station. 


In reversed flow regions, the streamwise convection terms in the x-momentum and 
energy equations can be evaluated using flow-dependent difference operators. These flow- 
dependent operators must model the proper zone of influence at every grid point. For 
example, a three-point backward difference operator may be used when u > 0, while a 
three-point forward difference operator would be used in the reversed flow regions. This 
has been avoided in the past because of the requirement of repeatedly sweeping through 
a reversed flow region. A popular approximation (called the FLARE 26 approximation) 
avoids the need for multiple sweeps through separation regions by setting the streamwise 
convection terms in the momentum and energy equations to zero in separated regions. 
These one-sweep solutions are achieved at the expense of assuming that the upstream 
influence in reversed flow regions is small. However, if the viscous algorithm is part of 
a viscous-inviscid interaction code, repeated sweeping of the viscous flow is required as 
part of the viscous-inviscid iteration. In this case, upstream influence in the reversed 
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flow region can be obtained during the repeated sweeping required by the viscous-inviscid 
iteration for little additional computational cost. 

Near and in reversed-flow regions, the boundary-layer equations axe solved in the 
inverse mode to avoid singular behavior at the separation point. While various researchers 
use 6*, here the wall shear stress r w and the wake centerline velocity u wc are used because 
they can be readily and efficiently implemented. The algorithm is modified to operate 
in the inverse mode by replacing the pressure term in the x-momentum equation with am 
expression containing r w or u wc . These relations are obtained by applying the x-momentum 
equation at the wall and wake centerline. For example, at a wall the x-momentum equation 
yields 


PUvt = faunVy )nVv L (3) 

The right-hand side of Eq. (3) is evaluated using first-order finite-difference approximations 
which allow the elimination of the pressure term ( PizSfp ) from the x-momentum equation. 
However, in the inverse mode the value of u at k = 2 , i.e., u 2 , appears in the difference 
equation at each and every k index. If this term is treated implicitly, then the matrix is 
tridiagonal with an additional column of nonzero coefficients 


/ 62 c 2 

J h &3 c 3 

f * <*4 ^4 C4 


\ 


u = d 




tokmax 


c kmax—l I 
bkmax / 


( 4 ) 


This augmented scalar tridiagonal matrix system is efficiently solved using an algorithm 
developed specifically for this application. The algorithm uses an LU decomposition to 
solve the augmented scalar tridiagonal system 24-35 and costs only 40% more than the 
standard scalar tridiagonal inversion. This approach results in an implicit treatment of 
the inverse boundary conditions. 


B. Unsteady Three-Dimensional Formulation 

As will be shown in Part II of this paper, the steady two-dimensional formulation just 
described has been quite useful. However, the extension of the space-marching boundary- 
layer approach to three dimensions is not straightforward because there is not a clearly 
defined marching direction. Our experience in three dimensions is that a time-relaxation 
approach to solving the boundary-layer equations is preferable. 27-28 

By relaxing the boundary-layer equations in time, flow-dependent difference opera- 
tors which automatically adjust to the flow direction can replace problem-dependent space 
sweeps. Hence, one relatively simple code can be easily applied to a wide variety of 
three-dimensional flows without any changes to the algorithm. An additional benefit of 
relaxation schemes is that they are also more robust and forgiving of approximate bound- 
ary conditions than space-marching schemes. With a space-marching scheme, all variables 
must be specified at inflow and some side boundaries, whereas simple outflow, reflection, 
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periodic, or zero-gradient boundary conditions can often be used with & relaxation numer- 
ical algorithm. For example, it would be very difficult to specify the flow quantities near 
a wing tip, as may be required by a space-marching scheme. However, with a relaxation 
algorithm the tip values can be brought up during the iteration process using, for exam- 
ple, a zero-gradient boundary condition. On the other hand, compared to space-marching 
algorithms, relaxation schemes can be relatively computationally expensive when used to 
simulate simple steady flows. This drawback can be minimized if care is taken to design 
the relaxation algorithm so that a steady state can be reached in few iterations. Also, since 
a relaxation boundary-layer scheme can be approximately half as expensive as a second- 
order space-marching scheme on a per-iteration basis, when included in a viscous-inviscid 
interaction algorithm the relaxation boundary-layer scheme should result in a more efficient 
code. Finally, the same code can also be used to simulate unsteady flows. 

1. Equation Set 

Neglecting surface curvature, the compressible boundary-layer equations for the un- 
steady, three-dimensional flow of a perfect gas can be written in £(x,z), i](x y y, z), £(x,z) 
coordinates as 

x-momentum equation 

p(u t + Uu( + Vu v + Wu() = —(3(pt$ s +P(C») + ( 5a ) 

z-momentum equation 

p( w t + Uw ( + Vw n + Ww() = ~/3(p(Cz + P(Cz ) + (p-w v rj y ) v rj y (5 b) 

perfect gas relation 

P = pT (5c) 

energy equation (H'= total enthalpy) 

p(B, + UH ( + VH, + WH() = ny + p, (5 d) 


continuity equation 


pt + (pu)t£ x + ( pu) v rj x + (pu) ( Cx + ( pv)nVy + (pw)({z + (pw)r,Vz + (pu>)<(z = 0 (5e) 


where U, V, and W are unsealed contravariant velocities 


’ U " 


■&+ 

txU+ 

tz to' 

V 

= 

Vt+ 

1JzU+ 

*7yV+ IJzW 

w 


m (t+ 

(zU+ 
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In the interest of simplicity, the boundary-layer equations and algorithm are presented 
here neglecting surface curvature. To date, we have included the terms to treat body- 
conforming orthogonal curvilinear coordinates . 29-31 Eventually, a general treatment of 
surface curvature should be included . 32 

2. Solution algorithm 

Once again, because the pressure gradient forcing terms in the momentum equations 
are treated as given functions, the boundary-layer equations are weakly coupled and can 
be solved sequentially at each time step. As a result, a semi-implicit algorithm can be used 
at each time or iterative step, yet only scalar-like uncoupled equations are solved. Spa- 
tial derivatives in the momentum (and energy) equations are approximated with implicit 
second-order-accurate central-difference operators in the direction normal to the body ( 77 ), 
and flow-dependent second-order accurate operators in the other two directions (£ and £). 
That is, central differencing is used in 77 , and upwind differencing is used in £ and £. For 
example, when differencing a term such as u^. a backward difference is used if the coeffi- 
cient pU is positive and a forward difference is used if the coefficient pU is negative. In the 
continuity equation, the £ and £ derivatives are approximated with second-order-accurate 
central-difference operators, and trapezoidal-rule integration is used in the 77 direction to 
yield v. 

The time- accurate algorithm is presented below using conventional operators (defined 
in the nomenclature) and the second-order-accurate upwind operator 6 

a, - (4H) v, (i^C) * (^H) », (!^) 

where a is the convective coefficient (e.g., pU). Using a notation that space-time indices are 
not written unless changed (e.g., u = u!? fc) j,u n+1 = u?^ 1 , etc.), the semi-implicit scheme 
is described below. 

From the x-momentum equation, the velocity u at the new time step is updated from 


p(V t u n+1 + U8{U n+1 + V£ n ti n+1 + WS^u) = -/?(£*%> + C*%>) + *?y£>jO i 7 ?y£?« n+1 ) ( 6 a) 

where ( U = £ x u + ( z w, etc...) and where if U <0, will actually operate on explicit 
data. As noted below, care is taken to maintain stability of terms that axe explicitly 
three-point differenced. This system of equations is easily solved by inverting scalar- 
tridiagonal matrices in 77 . In a similar way, w is updated from the z-momentum equation 
(U = £*u n+1 -f Z x w,etc...) 


p(V t w n+1 +U8{W n+1 + V8 v w n+1 +W8(w) = -/3(U^(P+ CzS^p) + rj y S v (prj y Sr,w n+1 ) (6b) 
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and H is updated from the energy equation ( U = £ z u n+1 + £ z w n+1 ,etc...) 


p(VtH n+1 +U6tH n+1 + VS v H n+1 + WS ( H ) 

= +w’)” + ‘}] + v,p- 


(6c) 


Finally, obtain p = p/T and integrate the continuity equation for v using updated values 
of u,w, and p. 


Vv(P v ) n+1 = ~ 1+ 2 ?n + n* s v(P u ) + CxSi(pu) (6d) 

+ £*t>t(P w ) + Vz8t,(pw) + C zS((pw ) + V t p] n+1 

The coefficients p and Pr are then evaluated using a turbulence model, here either the 
Cebeci algebraic turbulence model 33 (generalized to three dimensions) or the Baldwin- 
Lomax model. 34 

The algorithm just presented is second-order accurate in space and first-order accurate 
in time. A simple corrector sequence can be used to achieve second-order time accuracy and 
provide additional stability of the explicitly differenced terms following the same technique 
employed in the two-dimensional space-marching algorithm described earlier. However, 
if only a steady-state solution is required, enhanced stability and faster convergence is 
obtained by implicitly including the diagonal contribution of any explicit operator. In this 
case the algorithm takes on a similarity to a SLOR scheme. Finally, in many steady-state 
cases an assumption of constant total enthalpy can replace the partial-differential energy 
equation. 

If no streamwise separation is present, these equations may be solved with the pressure 
specified (direct mode). If streamwise separation occurs, the inverse mode must be used 
and this is implemented as outlined for the two-dimensional case. It is interesting to note 
that it is possible to solve the three-dimensional boundary-layer equations in a mixed 
direct /inverse mode. For example, the x-momentum equation may be solved in the inverse 
mode, while the z-momentum equation can be solved in the direct mode, or vice-versa. 
The authors have observed that the saddle-point behavior at the streamwise separation 
line may be avoided if just the momentum equation in roughly the free stream direction is 
solved in the inverse mode (i.e., the other momentum equation can be solved in the direct 
mode). 

This algorithm readily lends itself to vector computer architecture by arranging the 
computation of the inversion matrix coefficients and the actual inversion in a vectorized 
loop. To do this, one spatial direction must be chosen along which to create the vec- 
tors. It is convenient to choose the direction along the body roughly normal to the free 
stream. Then at each time step, along each £ plane, all the coefficient matrices are com- 
puted, stacked together, and fed to a vectorized version of a scalar tridiagonal (or aug- 
mented tridiagonal) solver. Using this approach, the algorithm requires approximately 2 
ps / point /iteration on a CRAY-XMP processor. 
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3. Results 


Detailed experimental data have been obtained for the flow over a 6:1 prolate spheroid 
at angle of attack. 35-38 The «oo=45 m/sec transitional flow over the ellipsoid at 30° angle of 
attack was computed using the relaxation boundary-layer code. The experimental pressure 
and near-surface particle traces are presented in Figs. 2-3, respectively. The experimented 
particle traces were obtained by computing the streamlines of particles released on the 
body (and required to remain on the body) given the experimental wall skin-friction data. 


The flow was computed by solving the boundary-layer equations, including the non- 
unity geometrical metric coefficients outlined by Wang, 29 Cebeci, 30 and others, on a 45 
(axial direction) x 50 (normal to the body) x 25 (circumferential direction) mesh. A cubic 
spline representation of the experimental pressure distribution was specified as the forcing 
function. In Figs. 4a-c computed surface streamlines are presented for three conditions: 
fully laminar, transition at the cross-flow separation line, and fully turbulent. From these 
three figures it is clear that this flow is very sensitive to the laminar/turbulent state 
of the boundary layer. Specifically, in the fully laminar and transitional computations, 
the streamlines indicate two counterrotating streamwise vortices (indicated by the co- 
alescence/dispersion/coalescence streamline patterns in the meridional direction), while 
in the fully turbulent case, only one fairly small streamwise vortex is indicated by the 
one streamwise coalescence line. By comparison with the experimental particle traces 
it is apparent that the experimental flow is neither entirely laminar nor turbulent, and 
that’ the transitional result is in the closest agreement to the experimental data. The 
simple transition model used here has underpredicted the extent of turbulent flow, and 
as a result a slightly larger separated region is predicted than is observed experimentally. 
However, it is clear that these flows can be computed routinely with the present code and 
a refined transitional model could be developed using this code. Figure 5 is a close-up of 
the captured stagnation point. As mentioned earlier, with the present algorithm there is 
no special treatment of the stagnation region since there is no requirement to march away 
from the stagnation point. 
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Fig. 2 Measured pressure distribution 37 on a 6:1 prolate spheroid at 30° angle 
of attack, Re = 7.2 x 10®, and Moo = 0.13. 



Fig. 3 Measured near-surface particle paths for the flow over the prolate 
spheroid of Fig. 2. 



M 




Fig. 4 Computed near-surface particle paths for the flow over the prolate 
spheroid of Fig. 2: a) completely laminar; b) transition at the cross-flow 
separation line; c) completely turbulent. 
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Fig. 5 Closeup of captured stagnation point of the flow over the prolate 
spheroid of Fig. 2. 

The transonic flow over a NACA 0012 wing with a 20® sweep and an aspect ratio of 3 
has been studied by Lockman and Seegmiller 39 over a range of free stream Mach numbers 
and angles of attack. At Moo = 0.826, a = 2°, a strong shock runs across the wing (Fig. 
6) which contributes to the formation of complex streamwise and cross-flow separation 
regions. Unfortunately, the experimental data are not of sufflcient detail to supply the 
input for a boundary-layer computation. However, this flow has also been computed 
using a Navier-Stokes algorithm 40 (called TNS) which does supply detailed results. In 
general, the pressure distribution generated by the TNS code is in good agreement with the 
experimental data. Consequently, the pressure distribution generated by the TNS code can 
be used as input for a boundary-layer computation, and this was done by operating in the 
inverse mode and iterating to the TNS pressure distribution. It is important to note that 
if this iteration procedure were to converge to exactly the specified pressure distribution, 
saddle-point behavior would be encountered at any streamwise separation line. To avoid 
this, the relaxation parameter u> which appears in the expression to iteratively update the 
shear stress (r” +1 = r” + u>(p£ — Pj>Jvs)) * s ma de proportional to the absolute value of the 
computed wall shear stress. (A complete description of the scheme used to update the wall 
shear stress, Eq. (8) with p* = ptns * is given in §A.l of Part II of this report.) In this way, 
the specified pressure distribution is accurately matched, but the saddle-point behavior is 
avoided since the iteration procedure never reaches exactly the specified pressure at the 
separation line. The experimental oil-flow and computed near-surface particle paths are 
compared in Figs. 7-8. As shown, details of the flow near the tip, the “mushroom” 
separation region, and (farther inboard) the continued deflection of the flow through the 
shock are all captured. Overall, even for this complex case, the agreement between the 
experimental and computed boundary-layer flows is quite encouraging. 



Fig. 7 Experimental oil-flow patterns 39 Fig. 8 Computed near-surface parti- 
for the flow over the NACA 0012 cle paths for the flow over the N AC A 
wing of Fig. 6. 0012 wing of Fig. 6. 



PART II: INTERACTION SCHEMES 


Three approaches to the viscous-inviscid interaction problem are briefly reviewed. 
First, a conventional full-potential/boundary-layer transpiration boundary- condition match- 
ing procedure will be discussed. Although this approach is useful and fairly flexible, it is 
our experience that it is difficult to apply in some cases. For this reason alternative inter- 
action procedures are being explored. 28 ' 41-43 Two of these, a vorticity-interaction scheme 
and the Fortified Navier-Stokes (FNS) approach, are presented here. Many alternative 
interaction procedures jure being studied. However, we feel that the two presented here 
are promising and a discussion of them will outline some of the basic requirements for the 
next generation of interaction procedures. 

A. Full-Potential/Boundary-Layer Interaction 


The goal of a viscous-inviscid interaction algorithm is to obtain converged and com- 
patible viscous and inviscid flow field solutions. To achieve this the viscous-inviscid iter- 
ation scheme must allow each flow to be influenced by the other, yet remain stable and 
computationally efficient even when this interaction is strong. The “classical” interaction 
approach is to obtain an approximation to the inviscid flow, extract some information from 
the inviscid flow (e.g., p) to feed to the viscous algorithm, and finally extract information 
from the viscous algorithm (e.g., 6 *) and input it to another estimate of the inviscid flow. 
This cycle (Fig. 9) is continued until the inviscid and viscous solutions are converged 
and compatible. The interface schemes for transferring information between the inviscid 
and viscous algorithms are critical and still a subject of active research. In the next two 
sections, two-dimensional examples of the more “classical” interfaces are discussed. 


1. Inviscid to Viscous Interface 


If the flow remains attached, the boundary-layer equations can be solved with pressure 
specified (the direct mode), where the pressure is obtained from the inviscid solution. The 
need to switch to the inverse mode near regions of reversed flow complicates the inviscid 
to viscous interface in two ways. First, the interface must “sense” when the viscous flow is 
about to separate and then change from the direct to the inverse mode. Second, since the 
inverse forcing functions (r w and u wc ) cannot be supplied directly by the inviscid algorithm, 
interaction algorithms must be developed which supply these forcing functions. These 
interaction algorithms should supply inverse forcing function distributions which result in 
the inviscid and viscous flow algorithms rapidly converging to compatible solutions in as few 
viscous-inviscid iterations as possible. Since the converged solutions will not depend upon 
the form of the interaction algorithms, only the requirement of computational efficiency 
need influence the design of these interface algorithms. 
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Fig. 9 “Classical” interaction scheme. 

Monitoring the evolution of the wall shear stress is a simple way of determining if a 
boundary layer is about to separate. In a typical boundary-layer computation it is rea- 
sonable to begin the computation in the direct mode, switching to the inverse mode only 
if r w falls below a prescribed small positive value. Once beyond the reversed flow regions 
it is possible to switch back to the direct mode. However, for simplicity, in the present 
work the viscous flow computation is continued in the inverse mode even if the flow has 
reattached. 

Once in the inverse mode, inviscid to viscous interface algorithms must supply inverse 
forcing function distributions. A variety of schemes for updating the wall shear stress 
have been studied. The integral boundary-layer equation 44 indicates that skin friction is 
a function of the pressure gradient 


0 / 1 . 
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0 du t 
u e dx 


(2 + H-M e 2 ) + 


d0_ 

dx 


(7a) 


Likewise, evaluation of the s-momentum equation at the wall (Eq. (3)) suggests a relax- 
ation equation of the form 


+ m- tti T * \ / \ 

{*ir) (7i) 

where = r” or r” +1 and w is a relaxation parameter. It was found that schemes based 
upon this type of equation worked well for some separated flows, such as incompressible 
diffuser flow or transonic airfoil flow with weak shocks. However, strong shock cases gener- 
ate large pressure gradients which caused all of the schemes based on pressure gradient to 
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be unreliable. A scheme based simply upon the difference between the viscous and inviscid 
pressures was found to be more reliable and equally fast. The results presented in this 
work were obtained using the simple scheme 

= K + - p? +1 ) (8) 

In all cases, an over-relaxation factor of ten (i.e., u> = 10 in Eq. 8) was used to accelerate 
convergence. 

This scheme is based on an observation about the nature of boundary-layer flow. It 
was noted that at a given streamwise station an increase in the specified r w will result in an 
increased computed u e or a decreased pressure (Fig. 10). Hence, if the pressure computed 
by the viscous algorithm is greater than the desired inviscid pressure, the r w should be 
increased from the value used in the previous iteration. A similar argument holds in the 
situation where an increased viscous pressure is desired, in which case a decrease in t w 
is required. Carter 45 developed a similar scheme to update the displacement thickness 
specified to an inverse code. As mentioned previously, the only criterion these schemes 
must meet is that they produce matched inviscid and viscous solutions in as few iterations 
as possible. 
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Fig. 10 Relation between r w and p at a given streamwise station. 


2. Viscous to Inviscid Interface 

The classical method of introducing the influence of the viscous flow upon the inviscid 
flow is to add the displacement thickness distribution to the original body shape, and 
then compute the inviscid flow about this new body. The advantage of this method is its 
simplicity. However, there are two disadvantages to this approach. 
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If the displacement thickness approach is used, the inviscid grid must be regenerated 
after each viscous iteration. The body-conforming inviscid grid must be regenerated be- 
cause the body about which the inviscid flow is computed changes when the displacement 
thickness distribution changes, after each viscous iteration. The computational cost of 
regenerating the inviscid grid may be minimized by simply shearing the original inviscid 
grid. However, we have found that this can still be costly and may lead to unsuitable 
inviscid grids. 

The second disadvantage of the displacement thickness approach is the possibility of 
supercritical viscous-inviscid interaction. First, it should be explained that subcritical and 
supercritical interactions are distinguished by the response of the flow deflection angle 0 
(the ratio of the inviscid velocity components normal and parallel to the body) generated 
by the viscous flow to the streamwise pressure gradient. In an adverse pressure gradient a 
subcritical interaction results in an increased 0, while a supercritical interaction results in 
a reduced 0. (This nomenclature was first introduced by Crocco and Lees. 46 ) Because of 
these properties a supercritical interaction cannot result in a smooth transition to separa- 
tion, a process in which the pressure and boundary-layer thickness both increase. Rather, 
before a supercritical interaction can produce separated flow, a transition or jump back 
to a subcritical condition must occur. This required jump or discontinuity corresponds to 
a saddle-point singularity and is similar to the throat singularity of one-dimensional noz- 
zle flow. Hence, any interaction procedure which leads to supercritical behavior and may 
encounter this singularity should be avoided. Le Balleur 47 has shown that the classical 
displacement thickness interaction becomes supercritical when M e « 1.5 — ► 2.0, depending 
on the history of the turbulent boundary layer. This implies that the classical displacement 
thickness approach should not be used in strong transonic viscous-inviscid interactions. 

An alternative to adding the displacement thickness distribution to the original body 
shape is to impose a transpiration velocity at the body surface, as first suggested by 
Lighthill. 46 In the transpiration velocity approach, the original inviscid flow tangency con- 
dition at the body surface is modified to enforce a velocity component normal to the 
body. This normal or transpiration velocity component deflects the inviscid flow away 
from the body, and thus simulates the displacement of the inviscid flow by the viscous 
flow momentum defect. The advantages of the transpiration velocity approach are that 
the inviscid grid need not be regenerated after each viscous iteration, and the interac- 
tion always remains subcritical. 47 Hence, the transpiration velocity approach is the more 
computationally efficient and robust viscous to inviscid interface approach. 

The transpiration velocity distribution should result in an inviscid streamline coincid- 
ing with the equivalent body shape obtained when the displacement thickness distribution 
is added to the original body shape. Such a transpiration velocity distribution may be 
computed as 46 


r. = (9) 

pe OS 

where p e and u e are the inviscid values at the body (or wake centerline), and in this 
instance s is the distance along the body (or wake centerline). Equation (9) is obtained 
by integrating the difference between the inviscid and viscous continuity equations from 
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y = 0 — * 8 (i.e., by integrating across the boundary layer), and simplifying the obtained 
result using two assumptions. These are the thin boundary-layer assumption, and the 
assumption that the inviscid flow is essentially uniform from y = 0 — ► 8. In the flow fields 
of interest in this work, these are reasonable assumptions. The transpiration velocities are 
then built into the body and wake boundary conditions of the inviscid code. In the case 
of the TAIR full-potential code 49 used in this work, Vs is converted to perturbations of 
the inviscid contravariant velocities (U p and V p ) which are included in the inviscid residual 
computations along the body and wake centerline (see Refs. 24-25 for details). 

3. Wake Treatment and Curvature Effects 

In many cases it is necessary to account for the viscous flow curvature and the pressure 
jump which occurs across curved streamtubes. Therefore, the method of accounting for the 
pressure variation across the shear layers developed by Lock and Firmin 50 is incorporated 
into the present algorithm. Suffice it to say that these corrections add a substantial degree 
of complexity to the code, but are important for high angle of attack cases and airfoils 
with a high aft-loading. 51 The reader is referred to Refs. 25 and 50-51 for details. 

4. Complete Interaction Procedure 

All of the components of the viscous-inviscid interaction which have been presented 
are shown schematically in Fig. 11. In an inviscid- viscous iteration it is not necessary, 
or even desirable, to iterate each equation set to convergence before updating the other. 
Indeed, the ratio of inviscid to viscous iterations (rf/ v ) has a significant influence upon the 
efficiency of the computation. If r{/ v is too large, then a strong viscous-inviscid interaction 
will not be resolved and instability will result. For example, if during the iteration process 
a separation bubble develops under a shock, the shock will often rapidly weaken and move 
upstream. If the viscous flow is not updated often enough, then the separation bubble 
will not “follow” the shock upstream, but will generate a spurious reacceleration of the 
subsonic inviscid flow downstream of the shock. That is, if is too large, a lag develops 
between the inviscid and viscous solutions which may result in a breakdown in the viscous- 
inviscid interaction and numerical instability. In contrast, if the viscous flow is updated 
more often than is necessary to resolve the viscous-inviscid interaction (r^ v is too small), 
only additional computational cost is incurred. Therefore, it is prudent to recompute the 
viscous flow after each inviscid iteration (r</ w = 1) when a flow field is first simulated or if 
it is known that a very strong viscous-inviscid interaction occurs. 

5. Results 

McDevitt 52 compiled the experimental peak Mach number on the 18% circular-arc 
airfoil with a trailing edge splitter plate as a function of free stream Mach Number. Figure 
12 compares the results of the present interaction method to this experimental data. The 
comparison with experiment is quite encouraging. Considering the wind tunnel wall effects 
documented by Levy 53 the slight discrepancy at the lower Mach numbers is expected. The 
discrepancy at the higher Mach numbers may be due to wind tunnel wall effects, but could 
also be due to the assumption of isentropic inviscid flow. The isentropic flow shock jump 
deviates significantly from the Rankine-Hugoniot flow shock jump as the pre-shock Mach 
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Fig. 11 Complete interaction scheme. 

number approaches 1.3. The dramatic difference between the inviscid and viscous peak 
Mach numbers is indicative of the strong viscous-inviscid interaction being modeled. 

Because detailed experimental results are available, 54 the RAE 2822 has become a 
popular supercritical airfoil for the verification of airfoil codes. Therefore, the M = 
0.730, Re oo = 6.5xl0 6 flow about the RAE 2822 airfoil at a Ci = 0.803 was computed. 
Figure 13a compares the experimental C p distribution, the C p distribution computed by 
the present method, and the C p distribution computed by Mehta 55 using a thin-layer 
Navier-Stokes code. There is significant wind tunnel wall interference in this data. Hence, 
the experimentalists suggest that computations not be performed at the measured angle 
of attack (a = 3.19°), but rather at a corrected angle of attack of approximately a = 2.79° 
or at the measured C\ = 0.803. Mehta performed his computations at the suggested 
corrected angle of attack (and computed a Ci = 0.793), while the present computations 
were performed at a = 2.81° to match the measured Ci = 0.803. It appears that some 
discrepancy is introduced because of the first order accurate resolution of the supersonic 
region by the inviscid method used in the present work. In particular, the details of the 
pressure distribution in the supersonic regions axe not resolved. Also, the Navier-Stokes 
computation was able to resolve the shock profile more accurately than the present method. 
This discrepancy is due in part to the high degree of clustering in the Navier-Stokes grid 
in the region of the experimentally observed shock location. However, overall the three 
pressure distributions are in good agreement. 

The experimental and computed C/\ e distributions for the flow over the upper surface 
of the RAE 2822 are presented in Fig. 13b. The boundary layers were tripped at 3% 
percent of chord in both the experiment and the two computations. Considering the 
difficulties involved in measuring skin friction, the agreement between the computed and 
experimental Cf | e distributions is fairly good. Both of the computations predict incipient 
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o EXPERIMENTAL DATA (McDEVITT) 



Fig. 12 Calculated and experimental M pe a k versus Moo for the Re = 8 x 10 6 
flow about an 18% circular-arc airfoil with a trailing-edge splitter plate. 

shock-induced separation, and the present method predicts a small extent of trailing edge 
separation. 

Figures 13c-d compare the computed boundary-layer displacement thickness and shape 
factor distributions with those determined experimentally for the flow over the upper sur- 
face of the RAE 2822. Because it is difficult to determine the boundary-layer thickness 
from Navier-Stokes results, the Navier-Stokes integral boundary-layer parameters are not 
smooth. Nevertheless, the very rapid growth of the displacement thickness under the shock 
and near the trailing edge are accurately predicted. These results suggest that the eddy- 
viscosity turbulence models are adequate for computing attached and mildly separated 
airfoil-type flow fields. Also, both of the computed shape factor distributions (Fig. 13d) 
predict the basic experimental trends. 

Chen et al. 5# extended this technology and a Navier-Stokes formulation 57 to the porous 
airfoil application shown in Fig. 14. (Porous airfoil research, see Ref. 56 for a complete 
reference list, has suggested that adding a porous surface to a wing might be a feasible drag 
reduction approach.) Some representative results are indicated in Figs. 15a-b. A number 
of interesting characteristics about this flow and the two formulations were observed during 
the work. In particular, it was found that for moderate cases the Navier-Stokes and viscous- 
inviscid interaction codes both performed well, with the interaction code being significantly 
more cost effective. However, for severe cases (very strong shocks or high blowing rates) 
it was necessary to rely on the Navier-Stokes formulation. It seems that the isentropic 
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o EXPERIMENT 




Fig. 13 Calculated and experimental results for an RAE 2822 airfoil at = 
0.730, .Reoo = 6.50 x 10®, Ci = 0.803: a) C p distributions; b) Cf\ e distributions; c) 
8* distributions; d) H distributions. 
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Fig. 14 Porous airfoil in transonic flow. 



x/c x/c 


Fig. 15 Comparison of C p distributions for a NACA 0012 airfoil at Moo = 
0.8, a = 0.0% and Re — 4.09 x 10® : a) computed and experimental results for a 
solid airfoil; b) computed results for a porous airfoil (porous over region of 
0.3 < x < 0.5). (TLNS-Thin Layer Navier-Stokes, IBL-Interactive Boundary 
Layer) 
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assumption (used in the full-potential formulation) and the boundary-layer assumptions 
are being violated. 


This formulation has also been applied to a tri-element augmentor wing configura- 
tion by Flores and Van Dalsem. 58 Figures 16a-b indicate representative results from this 
work. The flow about this configuration has also been computed using a Navier-Stokes 
formulation. 59 Overall, the agreement between the Navier-Stokes and viscous-inviscid in- 
teraction results and the experimental data is only fair. Furthermore, it is not at all clear 
that the Navier-Stokes results are superior to the viscous-inviscid results. However, the 
viscous-inviscid interaction computation required at least an order of magnitude less com- 
puter time. From this work one might conclude that the viscous-inviscid interaction code 
gives about the same accuracy for much less cost, and therefore is a better engineering 
tool. On the other hand, the viscous-inviscid interaction code is more complex because of 
the logic required to track the six boundary-layers and three wakes. 


Overall, it is felt that these results indicate that a “conventional” full-potential /bound- 
ary-layer interaction scheme is a practical and useful engineering tool for flows of moderate 
complexity. However, as either the flow physics or the geometry become more complex, it 
becomes much more difficult to apply a “ conventional” viscous-inviscid interaction formu- 
lation. In this case one might consider using a Navier-Stokes formulation, or alternative 
viscous-inviscid interaction schemes. 
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Fig. 16 Example results of application of viscous-inviscid interaction (TAUG-V) 
to the augmentor wing geometry at Moo = 0.7, a = 1.05°, and Re = 12.6 x 10® : a) 
C p distributions; b) Viscous velocity vectors in the throat region. 
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B. Vector-Potential/Boundary-Layer Interaction 

An advantage of conventional viscous-to-inviscid interaction using either transpiration 
or an effective body displacement is that the inviscid flow can be computed using a grid 
that is relatively coarse compared to what is needed for the viscous grid. Moreover, the 
inviscid grid need not necessarily be body conforming (e.g., a Cartesian inviscid grid can 
be used) in the way that a high Reynolds number Navier-Stokes grid must be. 

The disadvantage of conventional viscous-to-inviscid interaction is that features such 
as wakes require special treatment, and matching to a highly rotational inviscid outer flow 
is not straightforward. Our own experience with the previously described algorithm is that 
too much logic (such as the wake curvature correction) must be programmed. A thin-layer 
Navier Stokes code, while costly, is a less complex code which is more readily adapted to 
new problems. 

Instead of imposing the effect of the viscous flow on the inviscid flow through a bound- 
ary condition, the viscous flow effect can be embedded into the inviscid flow by means of a 
forcing function. In particular, the vorticity of the viscous flow can be embedded, through 
a field forcing term, directly into the “inviscid” algorithm. Such a viscous-inviscid inter- 
action procedure has some advantages over the approach described in the previous section 
and is briefly discussed in the next section (others are pursuing related ideas; e.g., Halim 
and Hafez, 42 and Whitfield 43 ). 

1. Vorticitv Transfer Interaction 

The Euler equations in nonconservative form can be related to the stream function 
and other potential-like decompositions of the velocity field and the vorticity. In this 
form, the Euler equations can be solved very efficiently. For example, the steady inviscid 
governing equations in three dimensions can be written in terms of the dual-potential 
velocity decomposition defined by 


q = V^ + VxB; B = 

or 

/ u\ / <t > 9 + *l>y - Xz\ 

q = ( V = [ - V>* (10) 

\wj \<l>z + X* ~ #v / 

where <f> is the usual potential function and B is the vector potential. Because four func- 
tions are used to define three velocity components, a constraint needs to be imposed on 
the potential functions to uniquely specify the velocity field. Subject to the consistency 
constraint 

V • B = i?* + Xy + — 0 (11) 

the governing equations take on the symmetric form 60 
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continuity 


(p<f>z)z + ( P<f>y)y + (p<f>z)x — ~ \{^y ~ Xz)p\z ~ [(^* ~~ tke)p]y — [(>l* 

-*»)?]» (12a) 

vorticity 

V 2 B = -w 

( 124) 

Crocco equation 

qxw = Vh 0 — TVs 

(12c) 

entropy convection 

US Z + V3 y + ws z = 0 

stagnation enthalpy convection 

(12 i) 

w(ho)« + v(ho)y 4- w(ho) z = 0- 

(12c) 

nonisentropic Bernoulli 


M(fc)-K£)]*- WM 

(12/) 

The vorticity appears explicitly in Eq. 12b. Hence, it is straightforward to impose a 
vorticity field generated by any source (e.g., viscosity) on the velocity field. In place of an 
effective displacement-thickness boundary condition, viscous flow effects can be embedded 
into the dual-potential inviscid flow equations as sketched in Fig. 17. The dual-potential 


equations are solved throughout, while in viscous flow regions, the boundary-layer equa- 
tions are also solved to resolve the vortidty, entropy, and stagnation enthalpy fields. Within 
the viscous layer, the inviscid convection equations for entropy and stagnation enthalpy 
axe discarded, and the boundary-layer vorticity is fed directly into the right-hand side of 
Eq. (12b). Likewise the entropy needed in the Bernoulli equation (Eq. (12f)) can be 
obtained from the boundary-layer equations, or the boundary-layer equations can directly 
supply density. Outside of the viscous layer, Eqs. (12a-c) can be used to account for vor- 
ticity in the flow, and although they do not generate vorticity, they can adequately convect 
vorticity. Like a Navier-Stokes approach, essentially one governing set of equations (the 
dual-potential form of the Euler equations) can be used throughout the flow field. The 
boundary-layer equations are also solved on a subset of the computational domain, but 
only to generate the vorticity, entropy, and stagnation enthalpy gradients. In this way, 
the inviscid and viscous equations do not interact on a boundary alone, but by embedded 
forcing functions interact over the entire viscous portion of the computational domain. 

As this approach allows vorticity transfer into the outer inviscid flow, it would seem 
to be a more general viscous-inviscid interaction than an effective displacement thickness 
boundary condition. However, for this procedure to work the inviscid mesh must be 
sufficiently refined to retain the concentrated vorticity near the wall. Therefore, if this 
procedure is to be practical, it must be possible to efficiently solve the global inviscid 
equations on a grid which is refined in the normal direction in the vicinity of the wall. The 
Poisson-like equations which arise with the dual-potential velocity decomposition can be 
efficiently solved on such a refined grid. 
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Fig. 17 Sketch of computational domains in which the dual-potential equations 
are solved throughout on a Navier-Stokes-like grid and the boundary-layer 
equations are also solved near the wall to generate vorticity. 

2. Results 

Preliminary applications of this vorticity interaction procedure have been reported 
by Steger and Van Dalsem 41 to a two-dimensional inlet /diffuser flow, and by Rao et al. 61 
for three-dimensional incompressible flow over a trough. In both test cases vorticity was 
confined to the viscous layer and did not convect into the main flow. 

The entire two-dimensional inlet/diffuser flow (indicated in Fig. 18), including the 
flow exterior to the inlet, was solved. On the lower surface of the inlet the grid was 
refined, and the boundary-layer equations were also solved (Fig. 18). The two-dimensional 
boundary-layer scheme described previously was used to solve the compressible boundary- 
layer equations at a low Mach number. Near and in regions of reversed flow the boundary- 
layer equations were solved in the inverse mode. Computed velocity profiles obtained from 
the “inviscid” equations are shown in Fig. 19. Figures 20a-c show the profiles in more 
detail, and also show the velocity profiles computed from the boundary-layer equations 
using Moo = 0.2 and Re = 4.4 x 10 5 based on the inlet length. Although the boundary-layer 
equations are solved on a finer grid, the agreement is excellent except in the region x = 0.6 
to x = 0.8 where small errors are observed. These errors may be partially attributed to 
interpolation error. (Vorticity was interpolated only as a function of y. The x-variation of 
the grid was ignored.) A typical separated profile shown in Fig. 21 demonstrates that if 
the boundary-layer vorticity is resolved, the inviscid equations can give good results even 
when too coarse an inviscid grid is used. The no-slip boundary condition was not imposed 
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on the Poisson equation, and the departure from no slip at the wall is an indication of any 
error. In this case the maximum slip velocity reached 2% of the free stream at the first 
profile shown in Fig. 20c. 



Fig. 18 Sketch of zonal domains in which the stream-function equations are 
solved everywhere and the boundary-layer equations are superimposed only 
on the inner lower surface of the inlet/diffuser wall. (Computations performed 
with yi — l and y a = 1.07.) 



the inviscid stream function with vorticity prescribed from the boundary-layer 
solution. 
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Fig. 21 Good agreement for a typical separated profile in which the boundary- 
layer resolved vorticity is embedded into the inviscid flow equations with in- 
adequate grid resolution. 

Using this scheme, three-dimensional calculations have also been carried out by Rao 
et al. 61 For the three-dimensional trough defined by 

z(x,y) = —0.03 sech(4a; — 10) sech(4y — 6) (13) 

solutions are presented in Figs. 22-23 which were obtained at an upstream reference 
Reynolds number of 8000. Figure 22 is a carpet plot of the displacement thickness. The 
vorticity components were evaluated from the three-dimensional boundary-layer algorithm 
just described. Comparison with other numerical results as obtained from Davis et al. 62 
are presented in Fig. 23. Figure 23 is a plot of the variation of the x-component of shear 
stress in the streamwise direction at the plane of symmetry. The agreement is relatively 
good. Edwards 63 obtained a solution for this geometry using an interacting boundary- 
layer algorithm coupled with an inviscid small-disturbance analysis, while Davis et al. 62 
obtained a viscous flow solution using a vorticity /stream-function type formulation for the 
Navier-Stokes equations. The present method worked well for attached flows, however 
some difficulty in obtaining converged solutions in separated flow regions was reported. 61 
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Fig. 22 Carpet plot of the predicted displacment thickness distribution for the 
three-dimensional trough geometry. 



Fig. 23 Comparison of the z-component of the shear stress in the streamwise 
direction along the trough centerline. 
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C. Fortified Navier-Stokes Scheme (FNS) 

As an alternative to the more familiar “zonal” concept of solving the various flow 
zones on separate grids and patching the zones together, the Navier-Stokes equations can 
be applied throughout and the solution of simplified or subset equations can be used to force 
or “fortify” the general or global algorithm (Fig. 24). By adding subset equation forcing 
terms to the Navier-Stokes equations, the Navier-Stokes and subset equations interact 
strongly over entire regions rather than just at interface boundaries. Moreover, the subset 
equations can be applied selectively to only those regions where they are clearly valid. 
Specifically, a boundary-layer solver (capable of quickly resolving thin shear layers) can be 
used to force a global Navier-Stokes scheme. However, if during the iteration procedure 
the boundary-layer approximations become suspect in a certain region, the boundary-layer 
to Navier-Stokes forcing can be turned off, and the region can be resolved with only the 
global Navier-Stokes scheme. Because of this flexibility, the generality of the Navier-Stokes 
equations is retained, while some of the efficiency of the subset algorithms is recovered. A 
method for applying this idea, the FNS approach, is described in this section. 




Fig. 24 Boundary-layer forced Navier-Stokes schemes: a) Zonal scheme in 
which the boundary-layer and Navier-Stokes solutions are matched at a com- 
mon boundary of a patched grid; b) Spatial-forcing scheme in which the 
boundary-layer equations are solved on an independent grid, and the bound- 
ary-layer solution is forced into the Navier-Stokes algorithm over the entire 
viscous region. 
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1. Concept Development 

To illustrate the forcing concept which will be used with the Navier-Stokes equations, 
it is useful to begin with a linear system of partial-differential equations 

d t Q + £Q=0 (14) 

where Q is a vector, £ is a linear spatial differential operator, and /? may contain terms 
from linearizing a nonlinear system of partial-differential equations. With finite- difference 
approximation the original nonlinear system of partial-differential equations can be re- 
placed by a large system of linear ordinary differential equations 

d t Q + LQ = b (15) 

where Q is now the vector of the variables at each grid point, L is a matrix, and b contains 
ft and boundary conditions. In general, Eq. (15) is then solved using any one of a variety 
of time or relaxation schemes. However, in the case of the Navier-Stokes equations, Eq. 
(15) can be a very large and stiff system of equations. 

Assuming that a solution Q/ is known in some region (for example, from the solution 
of a simpler but still valid set of equations) then this solution can be used to reduce the 
computational domain of the original problem. Alternatively, Q/ can be embedded into 
the system of equations by means of a forcing function such as 

d t Q + LQ = b + x(Qf - Q) (16) 

where the parameter x is set to a large positive value if Q / is known and is set to zero 
everywhere else. 

The steady-state solution of Eq. (15) is 

Q = L -1 b 

while the steady-state solution of Eq. (16) is 

Q = (L + xI)-‘(b + xQ/) 

For sufficiently large values of x, Q — *• Q /, while for x = 0 everywhere Q = L -1 b, and 
some blend of these solutions is obtained for moderate or variable values of x* 

For a constant value of x the exact solution of Eq. (16) is 

Q = ce-< l+ *»< + (L + *I)-‘(b + xQf) 

so that a time-accurate path to a steady state will be damped as e~^ L+xI ^. Hence as x is 
increased, the steady state is reached faster. 

For an iterative scheme in which at least xQ is treated implicitly, such as 

(1 + A X ) Q” +1 = (I - AL)Q" + A(b + XQ/) (17) 
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the iterative solution is 


q-m . (l^a ) " Q > + k V P-M-r (b + q v 

q U + fcJ q £'(l + kx) m+l{ XQ,) 


and if 


/ I — hL \ 

Vl + hX/ 


< 1 then 


(18) 


(L + xI)-> 


= k E 


(I - kL) m 

(1 + kx ) m+1 


This simple analysis indicates that both the transient solution decays more quickly 
and the series approximation to the steady-state solution converges more quickly as x is 
increased. Hence, if a solution Q/ is known, it can be used to quickly drive an iteration 
scheme to a steady state or damp errors in transient solutions. Finally, although x is de- 
noted as a parameter, it could be any positive definite operator that enhances convergence. 
The reader is referred to Ref. 64 for further details. 


2. An Application 

In a typical high- Reynolds number Navier-Stokes simulation the fine-grid resolution 
is generally provided in a thin zone near the body surface, and the outer flow is effectively 
resolved as rotational inviscid flow (unless the turbulent coefficients are quite large, in which 
case extensive modeling is required). Because of the fine-grid resolution required near the 
body, a given algorithm often operates much less efficiently on the Navier-Stokes equations 
than it does on the Euler equations, even though the work per step may be similar and 
the viscous terms may enhance stability. However, on this same refined viscous grid, the 
boundary-layer equations can be efficiently solved. As a result one can speculate that by 
using the FNS approach and a boundary-layer algorithm to supply Q/ in the shear layers, 
it may be possible to significantly improve the productivity of a Navier-Stokes algorithm. 

In this application, an implicit approximately factored algorithm with central differ- 
encing in the rj and £ directions and upwinding in the £ direction is used to solve the 
three-dimensional thin-layer Navier-Stokes equations. The basic two-factor solution algo- 
rithm is written as 65 


[/ + />«£(!+)" + hS ( C n - hRe- , S ( J- 1 M n J - D,| c l 

x + + AQ” = 

- A({«|[(F + )” - f + ] + [(§-)* - F- ) + MG” - Geo) + «<(H n - Hee) 

-■Re-‘I C (§” - See)} - (Z>.|, + U,|<)(Q” - Qoo) 


(19) 


Here h = At (first order in time), or fc = y , (second order in time) and a free stream 
base solution is subtracted out to improved accuracy in the far field. The operators 8 £ 

and 6^ are backward and forward three-point difference operators. The flux F has been 


35 



eigensplit and the matrices A y B,C , and M result from local linearization of the fluxes 
about the previous time level. Because central-space-difference operators are used in rj 
and £, implicit Di and explicit D e numerical dissipation terms are included in Eq. (19). 

There are a number of possible ways to add the forcing term to Eq. (19). The most 
straightforward way would be to add the term hx to either the first or second left-hand 
side factor, and the term hx(Qf — Q n ) to the right-hand side. However, this would lead 
to a factorization error on the left-hand side of the order of x^ 2 , which can become large 
as x >■ 1* A factorization which results in an error of the order of x -1 h 2 is given by 

[j(l + h X ) + hS b (A + ) n + hS ( d n - hRe- 1 6 < J~ 1 M n J - Di\< 
x [1(1 + hx)] -1 x [l( 1 + h X ) + hS{(A~) n + hS v B n - DiU] AQ n = (20) 

- At{*|[(F+) n - F+ ] + */[(F") n - F" ] + S v (G n - Goo) + $<(H n - Hoc) 

- iZe' 1 *<(§" - Soo)} - (D e \ v + D t | C )(Q" - Qoo) + MQ/ " Q n ) 

With this implementation, x can become very large without concern for large factorization 
errors. Another advantage of this form is that as x becomes large it contributes to the 
diagonal dominance of both left-hand factors. 

3. Results 

The FNS method was tested on a geometry which roughly simulates the separated 
flow on a swept infinite wing, as studied by Van den Berg and Elsenaar. 66 The flow was 
computed with both the standard thin-layer Navier-Stokes 65 and the FNS algorithms. The 
flow was first computed using the standard Navier-Stokes algorithm alone on both a fine 
mesh (29 points in the flow direction, 50 points normal to the wall, 5 points in the span 
direction) and a coarse mesh which has only 20 points in the critical normal direction. 
The same minimum normal spacing at the wall was used in both the coarse- and fine- 
mesh computations. This flow was also computed with the FNS approach using the coarse 
Navier-Stokes mesh (20 points in the normal direction) and, near the wall, a superimposed 
fine boundary-layer mesh (50 points in the normal direction). In all the FNS computations 
presented here, x was proportional to the vorticity; hence, it is automatically large within 
the boundary-layer and rapidly drops to zero near the edge of the boundary layer. The drag 
history in Fig. 25 shows that the coarse-grid standard Navier-Stokes computation does not 
predict the drag accurately, and that the FNS method obtains essentially the same drag 
level in 50 iterations that the fine-grid standard Navier-Stokes computation reached in 400 
iterations. The computed near-surface particle traces are shown in Figs. 26a-c. Both the 
FNS method and the fine-grid standard Navier-Stokes computations predict a constant 
chord-line separation line, whereas the coarse-grid standard Navier-Stokes computation 
does not quite capture this qualitative feature. Thus, a coarse- mesh standard Navier- 
Stokes computation should not be relied on to predict even the basic flow features. A 
detailed comparison of the FNS and fine-grid standard Navier-Stokes computations near 
the separation line shows that the FNS computation has predicted slightly larger turning 
angles. It is thought that this slight discrepancy is due at least in part to differences in 
the amount of numerical dissipation in the Navier-Stokes and boundary-layer algorithms. 
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ITERATIONS 

Fig. 25 Drag coefficient history versus iterations for the turbulent flow over 
a 7.75% sine wave bump with a 35° leading-edge sweep, Re = 5 x 10 5 , and 
Moo = 0.5. 



(c) 

Fig. 26 Computed near-surface particle traces for the flow of Fig. 25: a) 
Fortified Navier-Stokes; b) standard Navier-Stokes (fine mesh); c) standard 
Navier-Stokes (coarse mesh). 
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The swept infinite wing geometry described above was modified to yield a truly three- 
dimensional flow by reducing the aspect ratio of the wing to one. Also, to resolve spanwise 
changes, the grid dimension in this direction was increased from 5 to 15. The resulting 
particle traces (Figs. 27a-d) show the same trends as described for the infinite span ex- 
ample. Also shown in Fig. 27d is the result obtained when the viscous terms and no-slip 
boundary conditions are removed from the global numerical algorithm (thus making it an 
Euler formulation). In this case, the entire influence of viscosity must be carried by the 
boundary-layer algorithm, which is not a difficulty for this case. The drag history versus 
CRAY-XMP CPU time for these computations is presented in Fig. 28. As before, the 
coarse-grid standard Navier-Stokes computation is not accurate, and the fine-grid stan- 
dard Navier-Stokes computation is expensive, whereas the FNS (and Eider) computations 
yield the same drag level as does the fine-grid standard Navier-Stokes computation, but 
for one-twentieth of the cost. 






Fig. 27 Computed near-surface particle traces for the turbulent flow over a 
7.75% sine wave bump with a 35° leading-edge sweep, Re = 5 x 10 s , = 0.5, 

and AR=1: a) Fortified Navier-Stokes; b) standard Navier-Stokes (fine mesh); 
c) standard Navier-Stokes (coarse mesh); d) Fortified Euler. 
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CRAY— XMP CPU mc 

Fig. 28 Drag coefficient history versus CRAY-XMP processor time for flow of 
Fig. 27. 

CONCLUSIONS 

A variety of viscous-inviscid interaction schemes have been presented. A “classical” 
full-potential/boundary-layer interaction has been applied to the transonic separated flow 
about solid and porous airfoils, and a tri-element augmentor wing. This scheme is rep- 
resentative of the majority of viscous-inviscid interaction codes presently in use. These 
codes have proven very useful for moderately complex, mildly separated, two-dimensional 
flows (e.g., the augmentor wing application), and fairly simple three-dimensional flows. 
However, it is the authors* opinions that it will be difficult to extend the effective displace- 
ment thickness class of interaction schemes much beyond present capabilities. The normal 
pressure gradient and free shear-layer tracking logic needed by these schemes requires a 
great deal of attention for each new problem and is quite tedious for truly complex flows. 

The vorticity interaction scheme presented here is one of a number recently developed 
specifically to overcome the deficiencies of the effective displacement thickness interaction 
schemes. In this case, the displacement thickness concept is not used. The outer “inviscid” 
algorithm is tightly coupled to the viscous flow over the entire viscous flow domain. If 
necessary, the inviscid flow can convect vorticity downstream (i.e., as in a wake). To 
date, this approach has been applied to separated, but fairly simple, flows. However, it is 
straightforward to include the effect of normal pressure gradient in this scheme, and free 
shear-layer tracking is automatic. Hence some of the deficiencies of the “classical” effective 
displacement thickness interaction have been avoided. 
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The Fortified Navier-Stokes scheme is the next step toward complete generality. In this 
case, a general global algorithm is “fortified” (i.e., the accuracy and/or convergence rate are 
improved) in any region where a less complete, but more efficient, numerical formulation 
can be used. For example, the boundary-layer equations can be used to resolve thin shear 
layers. It has been shown that this approach can significantly improve the efficiency of the 
overall algorithm. 
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